/*
  Copyright (C) 2011, Kyungjoo Kim
  
  This file is part of LINAL (LINear ALgebra)
  
  All rights reserved.

  This library is free software; you can redistribute it and/or
  modify it under the terms of the GNU Lesser General Public
  License as published by the Free Software Foundation; either
  version 2.1 of the License, or (at your option) any later version.
  
  This library is distributed in the hope that it will be useful,
  but WITHOUT ANY WARRANTY; without even the implied warranty of
  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  Lesser General Public License for more details.
  
  You should have received a copy of the GNU Lesser General Public
  License along with this library; if not, write to the Free Software
  Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301  USA
  Also add information on how to contact you by electronic and paper mail.
  
  
  Kyungjoo Kim
  iamkyungjoo@gmail.com
*/
#ifndef LINAL_UTIL_FLOP_HXX
#define LINAL_UTIL_FLOP_HXX

#define LINAL_MUL(is_complex) ( (is_complex) ?  (6.0) : (1.0) )
#define LINAL_ADD(is_complex) ( (is_complex) ?  (2.0) : (1.0) )
#define LINAL_MEM(is_complex) ( (is_complex) ? (16.0) : (8.0) )

namespace linal {

  /*! FLOP counting - From LAPACK working note #41
   */
  inline Double_ get_flop_gemm( Int_ is_complex, Int_ mm, Int_ nn, Int_ kk ) {
    Double_ m = (Double_)mm;    Double_ n = (Double_)nn;    Double_ k = (Double_)kk;
    return ( LINAL_MUL(is_complex)*(m*n*k) +
             LINAL_ADD(is_complex)*(m*n*k) );
  }
  inline Double_ get_flop_syrk( Int_ is_complex, Int_ kk, Int_ nn ) {
    Double_ k = (Double_)kk;    Double_ n = (Double_)nn;
    return ( LINAL_MUL(is_complex)*(0.5*k*n*(n+1.0)) +
             LINAL_ADD(is_complex)*(0.5*k*n*(n+1.0)) );
  }
  inline Double_ get_flop_trsm_lower( Int_ is_complex, Int_ mm, Int_ nn ) {
    Double_ m = (Double_)mm;    Double_ n = (Double_)nn;
    return ( LINAL_MUL(is_complex)*(0.5*n*m*(m+1.0)) +
             LINAL_ADD(is_complex)*(0.5*n*m*(m-1.0)) );
  }
  inline Double_ get_flop_trsm_upper( Int_ is_complex, Int_ mm, Int_ nn ) {
    Double_ m = (Double_)mm;    Double_ n = (Double_)nn;
    return ( LINAL_MUL(is_complex)*(0.5*m*n*(n+1.0)) +
             LINAL_ADD(is_complex)*(0.5*m*n*(n-1.0)) );
  }
  inline Double_ get_flop_trmm_lower( Int_ is_complex, Int_ mm, Int_ nn ) {
    Double_ m = (Double_)mm;    Double_ n = (Double_)nn;
    return get_flop_trsm_lower( is_complex, m, n);
  }
  inline Double_ get_flop_trmm_upper( Int_ is_complex, Int_ mm, Int_ nn ) {
    Double_ m = (Double_)mm;    Double_ n = (Double_)nn;
    return get_flop_trsm_upper( is_complex, m, n);
  }
  inline Double_ get_flop_lu( Int_ is_complex, Int_ mm, Int_ nn ) {
    Double_ m = (Double_)mm;    Double_ n = (Double_)nn;
    if (m > n)
      return ( LINAL_MUL(is_complex)*(0.5*m*n*n-(1.0/6.0)*n*n*n+0.5*m*n-0.5*n*n+(2.0/3.0)*n) +
               LINAL_ADD(is_complex)*(0.5*m*n*n-(1.0/6.0)*n*n*n-0.5*m*n+        (1.0/6.0)*n) );
    else
      return ( LINAL_MUL(is_complex)*(0.5*n*m*m-(1.0/6.0)*m*m*m+0.5*n*m-0.5*m*m+(2.0/3.0)*m) +
               LINAL_ADD(is_complex)*(0.5*n*m*m-(1.0/6.0)*m*m*m-0.5*n*m+        (1.0/6.0)*m) );
  }
  inline Double_ get_flop_chol( Int_ is_complex, Int_ nn ) {
    Double_ n = (Double_)nn;
    return ( LINAL_MUL(is_complex)*((1.0/6.0)*n*n*n+0.5*n*n+(1.0/3.0)*n) +
             LINAL_ADD(is_complex)*((1.0/6.0)*n*n*n-        (1.0/6.0)*n) );

  }
  inline Double_ get_flop_ldl( Int_ is_complex, Int_ nn ) {
    Double_ n = (Double_)nn;
    return ( LINAL_MUL(is_complex)*((1.0/3.0)*n*n*n + (2.0/3.0)*n) +
             LINAL_ADD(is_complex)*((1.0/3.0)*n*n*n - (1.0/3.0)*n) );

  }
  inline Double_ get_flop_qr( Int_ is_complex, Int_ mm, Int_ nn ) {
    Double_ m = (Double_)mm;    Double_ n = (Double_)nn;
    if (m > n)
      return ( LINAL_MUL(is_complex)*(m*n*n-(1.0/3.0)*n*n*n+m*n+0.5*n*n+(23.0/6.0)*n) +
               LINAL_ADD(is_complex)*(m*n*n-(1.0/3.0)*n*n*n+   +0.5*n*n+( 5.0/6.0)*n) );
    else
      return ( LINAL_MUL(is_complex)*(n*m*m-(1.0/3.0)*m*m*m+2.0*n*m-0.5*m*m+(23.0/6.0)*m) +
               LINAL_ADD(is_complex)*(n*m*m-(1.0/3.0)*m*m*m+    n*m-0.5*m*m+( 5.0/6.0)*m) );
  }
  inline Double_ get_flop_q( Int_ is_complex, Int_ mm, Int_ nn, Int_ kk ) {
    Double_ m = (Double_)mm;    Double_ n = (Double_)nn;    Double_ k = (Double_)kk;
    return ( LINAL_MUL(is_complex)*(2.0*m*n*k-(m+n)*k*k+(2.0/3.0)*k*k*k+2.0*n*k-k*k-(5.0/3.0)*k) +
             LINAL_ADD(is_complex)*(2.0*m*n*k-(m+n)*k*k+(2.0/3.0)*k*k*k+    n*k-m*k+(1.0/3.0)*k) );
  }
  inline Double_ get_memory( Int_ is_complex, Int_ mm, Int_ nn ) {
    Double_ m = (Double_)mm;    Double_ n = (Double_)nn;
    return LINAL_MEM(is_complex)*(m*n);
  }
}
#endif
